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I. INTRODUCTION 

Dielectric cavities have recently attracted considerable 
attention due to the fabrication of microlasers jl], ^J. 
Various shapes have been studied both experimentally 
and theoretically: deformed spheres || ||, ||, deformed 
disks § | J, fij, 1 1 ill, squares |j§ and 
hexagons [ [14[ FTa. An efficient numerical strategy to 
compute optical properties of effectively two-dimensional 
dielectric cavities with more complex geometries is the 
subject of the present paper. 

Maxwell's equations simplify to a two-dimensional (re- 
duced) wave equation |16| 



VV = n 2 (r)fc 2 V> 



(1) 



with coordinates r = (x,y) — (r cos#,rsin#), piece- wise 
constant index of refraction n(r), (vacuum) wave num- 
ber k = lo/c, angular frequency u> and speed of light in 
vacuum c. In the case of TM polarization, the complex- 
valued wave function ip represents the z-component of 
the electric field vector E z (r,t) = Re[V>(r) exp (— iwt)] 
with i 2 = —1, whereas for TE polarization, ip represents 
the z-component of the magnetic field vector H z . 

The boundary conditions at infinity are determined by 
the experimental situation. In a scattering experiment 
the wave function is composed of an incoming plane wave 
with wave vector k and an outgoing scattered wave. The 
wave function has the asymptotic form (in 2D) 

V ~ V>in + Vw = exp (ikr) + f(6, k) 6XP { ' kr) , (2) 

where k = |k| and f(0, k) is the angle-dependent differ- 
ential amplitude for elastic scattering. In lasers, however, 
the radiation is generated within the cavity without in- 
coming wave, 



1p ~ Ipo 



h(0,k) 



exp (ikr) 



(3) 
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This situation can be modelled by a dielectric cavity with 
complex- valued n leading to steady-state solutions of the 
wave equation ([!]) . Alternatively, one can use real- valued 
n leading to states that are exponentially decaying in 
time. The lifetime r of these so-called resonant states or 
short resonances is given by the imaginary part of the 
wave number as r = — l/2clm(fc) with Im(fc) < 0. r is 
related to the quality factor Q = Re(w)r. The resonant 
states are connected to the peak structure in scattering 
spectra; see jl7| for an introduction. Resonant states 
have been introduced by Gamow p8| and by Kapur and 
Peirles @. 

The wave equation (Q) with the outgoing-wave condi- 
tion (^) can be solved analytically by means of separation 
of variables only for special geometries, like the isolated 
circular cavity (see e.g. Ref. ^0|) and the symmetric 
annular cavity |12j. In general, numerical methods are 
needed. Frequently used are wave-matching methods 0. 
The wave function is usually expanded in integer Bessel 
functions inside the cavity and in Hankel functions of first 
kind outside, so that the outgoing-wave condition (^) is 
fulfilled automatically. The Rayleigh hypothesis asserts 
that such an expansion is always possible. However, it 
can fail for geometries which are not sufficiently small 
deformations of a circular cavity pT| . It should be men- 
tioned that for a different kind of boundary conditions 
at infinity, the wave-matching method can work well for 
special strongly noncircular geometries, e.g. rectangular 
integrated microresonators [22 ], 

More flexible are, for example, finite-difference meth- 
ods; see e.g. Q. These methods involve a discretization 
of the two-dimensional space, which is a heavy numeri- 
cal task for highly-excited states. An even more severe 
restriction is that it is impossible to discretize to infin- 
ity. One has to select a cut-off at some arbitrary distance 
from the cavities and implement there the outgoing- wave 
condition (^). For these reasons, finite-difference meth- 
ods are not suitable for computing resonances in dielec- 
tric cavities. 

A class of flexible methods with better numerical effi- 
ciency are boundary element methods (BEMs). The cen- 
tral idea is to replace two-dimensional differential equa- 
tions such as Eq. (Q) by one-dimensional boundary inte- 
gral equations (BIEs) and then to discretize the bound- 
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aries. BEMs have been widely applied to geometries 
with Dirichlct boundary conditions (wave function van- 
ishes), Neumann boundary conditions (normal deriva- 
tive vanishes) and combinations of them [^3[ ^j, |I| M . 
Bounded states have been calculated in the context of 
quantum chaos; for an introduction see Refs. [|27| 28 1. 
For scattering problems consider, for example, Ref . [ 29 1 . 
Resonances have been computed for scattering at three 
disks by Gaspard and Rice p0| |. 

The boundary conditions for dielectric cavities, how- 
ever, are of a different kind: the wave function and its 
(weighted) normal derivative are continuous across a cav- 
ity boundary. An analogous quantum problem in semi- 
conductor nanostructures has been treated by Knipp and 
Reineckc |3l| . Their BEM is for bounded and scattering 
states only. The aim of the present paper is to extend 
their approach to resonances in dielectric cavities for TM 
and TE polarization, including a discussion of spurious 
solutions, treatment of cavities with symmetries and cav- 
ities with corners. 

The paper is organized as follows. The BIEs are de- 

ue 

Section III describes the BEM. Section Fvl 



rived in the framew ork of the Green's function technic 
in Sec. £| 



demonstrates the efficiency of the method with an exam- 
ple of two coupled hexagonal resonators. Finally, Sec. |v| 
contains a summary. 



II. BOUNDARY INTEGRAL EQUATIONS 

In this section we derive the BIEs for the general case 
of J — 1 optical cavities in an outer unbounded medium. 
As illustrated in Fig. [j], the space is divided into J re- 
gions Clj, j = 1, 2, . . . , J, in each of which the index of 
refraction n(r) = rij is uniform. Without loss of general- 
ity nj is set to unity, i.e. the environment is vacuum or 
air. We first concentrate on TM polarization where both 
the wave function ip and its normal derivative are contin- 
uous across an interface separating two different regions. 



To reduce the two-dimensional differential equation ([j]) 
to one-dimensional integral equations, we first introduce 
the Green's function, which is defined as solution of 



[V 2 + n 2 /fc 2 ]G(r, r'; k) = S(r - r') 



(4) 



where 8(r — r') is the two-dimensional Dirac 5-function, 
r and r' are arbitrary points within Qj. The outgoing 
solution for the Green's function is 



G{vy-k)^- l l H^\n k\v-v'\) 



(5) 



is the zeroth order Hankel function of first kind |?2] . 
Multiplying the ^-equation (|l|) by G(r,r'; k) and sub- 
tracting the resulting equation from the G-equation (^) 
multiplied by ip{ r ) gives 

^(r)J(r-r') = ^>(r)V 2 G(r, r'; k) - G(r, r'; fc)V 2 V(r) 

= V[V>(r)VG(r,r';fc)-G(r,r';fc)VV>(r)] . 

Integrating this equation over the region £lj yields on the 
l.h.s. ip(r') since r' S fij. Applying Green's theorem, the 
integral on the r.h.s. can be expressed by a line integral 
along the boundary curve Tj — dflj , such that 

V>(r')= I ds[i/;(s)d 1/ G(s, r'; k) - G(s, r'; k)d u tp(s)] . 

Note that the boundary curve may consist of a number of 
disconnected components T, — r^UT^U. . . as depicted 
Each component is assumed to be oriented 
counterclockwise, smooth, and not to be a part of flj 
itself, i.e. £lj is an open set. d v is the normal derivative 
defined as d v — ^(r)V| r ; v(v) is the outward normal unit 
vector to Tj at point r; s = s(r) is the arc length along 
Tj at r. The derivative of the Green's function is given 
by 



where h[ is the first order Hankel function of first 
kind fil and 



cosaiif Vjfclr-r'l) , (7) 




FIG. 1: Geometry and notation for the BIEs. The cavity with 
domain Q± is bounded by the curve Ti, the one with domain Q2 
is bounded by T2- The domain C3 is "bounded" by = Ti, 



-.(2) 



T2 and by a circle Too at a large distance. 



COS a 



"Ml 



(8) 



The limit r' — > Tj in Eq. (||) is not trivial since both the 
Green's function and its normal derivative are singular at 
r' = r. However, it can be shown that these singularities 
are integrable for smooth boundaries. This is obvious for 
the second part of the integral kernel in Eq. (||) since for 
small arguments z = rij fc|r — r'| 



HW(z)~-]nz 

7T 



(9) 



The first part is also integrable. At first glance, this 
seems to be surprising because for small arguments 



2i 

7T2 



(10) 
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However, this singularity is compensated by 



cos a 



1 , 



(11) 



where k is the curvature of the curve Tj at r(s), which 
is finite for a smooth boundary. The limit r' — > in 
Eq. (g) can be performed in the sense of Cauchy's prin- 
cipal value, see e.g. Ref. |27| , giving 

~V(r') = P / da[^(a)ft,G(«, r'; k) - G(s, r'; fc)0^(s)] . 

^ (12) 
Comparing the l.h.s of Eqs. (p) and (PL2) shows that r' S 

Tj gives the "average" of the results for r' S fij and 

r' G $1^ with i ^ j. 



For each region fij there is an equation as Eq. (12). 
Special attention has to be paid to the unbounded region 
Qj. It is convenient to consider instead a finite region 
bounded by a circle with a very large radius r as 
sketched in Fig. [j]. We distinguish three cases in the 
following subsections. 



A. Bounded quantum states 

The case of bounded states in the quantum analogue 
has been studied by Knipp and Reinecke pi] . Then, rijk 
has to be replaced by [2m(E — Vj)] 1 / 2 /h, where E is the 
energy, Vj with j = 1, . . . , J is a piece- wise constant po- 
tential, and h is Planck's constant divided by 2tt. The 
wave function and its normal derivative (weighted with 
the inverse of the effective mass m) are continuous at do- 
main boundaries. If E < Vj then the state is bounded, 
the wave function and its gradient vanish exponentially 
as r — > oo. Moreover, with Im(fc) = the Green's func- 
tion (||) vanishes as either r or r' goes to infinity. As a re- 
sult Too does not contribute to any of the BIEs. Note that 
Eq. (|l|) does not permit bounded states since n 2 k 2 > 0. 

Using the same notation as Knipp and Reinecke |3l[] 
we reformulate Eq. (|l2|) as a linear homogeneous BIE 



j ds[B(s',s)^(s) + C(s',s)tp(s)]=0 , (13) 

with B(s', s) = -2G(s, s'; k), C(s' ', s) = 2d„G{s, s'; k) - 
8{s — s'), and 4>(s) — d^ip(r). The entire set of BIEs can 
be written in a symbolic way as 



B 2 



c 2 



M 



o 



(14) 



V Bj Cj J 



where Bj and Cj represent the integral operators in re- 
gion rij . The lower half of the vector (</>, i/>)* contains the 
values of the wave function on the boundaries, and the 
upper half contains the values of the normal derivative. 
Note that each boundary curve has two contributions to 
Eq. jl4| ) with identical -0, (which are continuous across 
the boundary) but different Bj, Cj. 



B. Plane- wave scattering 

The scattering states in the related quantum problem 
have been discussed again by Knipp and Reinecke |sT| . 
In contrast to the case of bounded states, their results 
also apply to dielectric cavities. 

In region £lj the wave function has the asymptotic 
form as in Eq. (Q). The incoming wave i[>i n satisfies 
Eq. (0). Thus, ip can be replaced by ip — tpi n in Eq. M) 
giving 

0(r') = exp(ikr')+/ ds{[^(s) - xp in (s)]d u G(s, r'; k) 



(15) 



-G(s,r';fcM S )-0 in ( S )]} , 



where ipi n (s) = exp (ikr) and 4>i n {s) — ik^(r) exp (ikr) at 
r = r(s). The circle at infinity does not contribute to the 
BIE ( J15| ) as in the case of bounded states. The reason, 
however, is different as we shall see in greater detail in 
the following subsection when considering resonances. 

If r' is taken from the boundary then Eq. (|l5| ) can be 
written as inhomogeneous integral equation 



ds[B(s',s)^{s)+C(s',s)ij(s)} = 



ds[B(s',s)cl> in (s) + C{s',s)ij in (s)} 



(16) 



Together with the other J — 1 BIEs, which are of the 
same form as in Eq. (13), the resulting inhomogeneous 
system of equations is 



M 



1> 



M 



(17) 



with 



M a = 



( \ 



V Bj Cj ) 



(18) 



Having determined the solutions tp and <fi we can com- 
pute the differential scattering amplitude by evaluating 
Eq. ( |l5| ) for large r' and comparing the result with Eq. (^) 
giving 



1 + i f 

f( e > k ) = -—F=r f d« exp [-ik/r(s)] 
Ay/irk Jr , 



(19) 



Wnk 

{ikfV(s)[lp(s) - V'in(s)] + 4>{s) - 0in(s)} , 

where k/ = (k cos 9, k sin 9) and 9 is the detection angle. 
Here, \f(9, k)| 2 is the differential scattering cross section. 
The total cross section cr(k) = J d9\f(9, k)| 2 can be easily 
calculated from the forward-scattering amplitude, k/ = 
k = (k cos <j>, k sin (f>), with the help of the optical theorem 
(see, e.g., Ref. |17 ) 



<7(k) = 2 A /-Im[(l-i)/(0 = 6k)] 



(20) 
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C. Resonances 

We now turn to the BIEs for resonances. Comparing 
the scattering boundary condition (^|) and the outgoing- 
wave condition (|J) indicates that we possibly can use the 
scattering approach neglecting the incoming wave, that 
is Eq. ( |l7| ) with Mq — 0. Apart from the fact that k is 
now a complex number, this is then identical to Eq. Jl4| ) 
for bounded states. There is, however, one problem. The 
circle at infinity, Too, may give a nonvanishing contribu- 
tion 

100(0 = / ds[^(s)d v G(s, r'; k) - G(s, r'; k)d^(s)} 

(21) 

to the r.h.s. of Eq. M) because with lm(k) < neither 
the wave function ([|) nor the Green's function (||) van- 
ish at infinity. Gaspard and Rice |lo| have shown for a 
Dirichlet scattering problem that nonetheless loo(r') = 
if r' is at one of the scatterers' boundaries or if r' is 
at a large distance from these boundaries. We have to 
extend their result because (i) the problem of dielectric 
cavities involves a different kind of boundary conditions; 
(ii) we are interested in the wave function ip( r ') also m 
the near-field. We start with recalling that the circle at 
infinity, Too, is defined by r = const with r — > 00. Using 
the asymptotical behaviour of Hankel functions of first 
kind Ip3| 



exp[i(,z-m7r/2- tt/4)] (22) 



as z = k\r — r'\ — > 00, it can be shown that the Green's 
function in Eq. (^|) is asymptotically given by 

G(r,v';k)^g(d-9'yf-^l, (23) 



with 

g(d-e',r') 



1 + i 



exp \-ikr' cos (9 - 9')] . (24) 



Equation (g3|) has the same r-dependence as the 
outg oing - wave condition (||) . With G and if) appearing in 
Eq. (21) in an antisymmetric way it follows Joo(r') = 
for all r' £ ftj U Tj. The fact that /oo(r') vanishes for 
r G 1 j means that the BIEs (0) can indeed be used 
to determine the resonant wave numbers k. Moreover, 
since Ioo( r ') = also for r' E Q,j Eq. (||) can be used to 
compute the corresponding wave functions in the entire 
domain. 

Having established that the resonances are solutions 
of the BIEs ( [l4"| ) with complex-valued k, we now demon- 
strate that the BIEs (|l4|) posses additional solutions 
which do not fulfil the outgoing- wave condition (||). We 
study this in an elementary way for a single cavity of ar- 
bitrary shape. Outside this cavity sufficiently far away 
from its boundary, a solution of wave equation (Til) can 



- a 



$H$(kr)]exp(im9) 



be expressed as 

00 

V(r,0)= £ [aWH<V(kr)- 

m=— 00 

(25) 

with Hankel functions of first and second kind |32| and 
with unknown complex- valued parameters a^m and 01$ . 
Without boundary conditions at infinity, solutions as in 
Eq. (E5T) exist for all values of k. Boundary conditions 

— (2) 
that fix all parameters a m give rise to a discrete spec- 
trum of k; for instance, the outgoing-wave condition (||) 
requires d£) = for all m. Inserting the expansion (I2T 
into Eq. flfll) leads to 



(26) 



/oo(r') = 2 J2 o£ , 4(*r')fflpM) 



m— — 00 



Hence, /oo(r') vanishes identically for all r' € f2j U Tj 

(2) 

only in the case of a resonance, where a m = for all m. 

However, the circle at infinity does not contribute to 
the BIEs ( |I4| ) already if the weaker condition lco( r ') = 
for r' £ Tj is satisfied. We insert this condition into the 
l.h.s. of Eq. (|2^) and note that the r.h.s. is an expansion 
of a solution of wave equation (|l|) inside the cavity with 
"wrong" index of refraction n = nj = 1. The result is 
that the BIEs ( [l4| ) possess undesired solutions, namely 
bounded states of an interior Dirichlet problem, in addi- 
tion to the resonances. As one consequence, the solutions 
of the scattering BIEs (|l7|) are not unique whenever k is 
a solution of the interior Dirichlet problem. Note that 
this nonuniqueness has not been discussed by Knipp and 
Reinecke plj . 

A related problem is known for cases with Dirichlet or 
Neumann conditions; see, e.g., Refs. |23, 25|. There have 
been several attempts to modify the BIEs in order to get 
rid of these "spurious solutions" . Some of these modi- 
fications could be applied to the present case, but this 
would result in singular integrals which are hard to deal 
with numerically. Fortunately, the spurious solutions are 
not a severe problem for our purpose. We can distinguish 
them, in principle, from the resonances in which we are 
interested in. The former have Im(fc) = whereas the 
latter have Im(fc) < 0. 



D. TE polarization 

In the case of TE polarization, Eq. (|l|) is valid with ip 
representing the magnetic field H z . The wave function 
-0 is continuous across the boundaries, but its normal 
derivative is not, in contrast to the case of TM polariza- 
tion. Instead, n(r)~ 2 d u ip is continuous [ fl6| . 

This new boundary condition can be easily incorpo- 
rated in the BEM by defining (j> = n~ 2 d,ji{j, B(s',s) — 
— 2G(s, s'; k)n 2 and <fii n accordingly in equations like 
Eqs. (|T^) and ([l6|). We remark that the spurious so- 
lutions are not affected by this change of boundary con- 
ditions. 
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E. Symmetry considerations 



III. BOUNDARY ELEMENT METHOD 



Many dielectric cavities studied in the literature pos- 
sess discrete symmetries. For example, the elliptical cav- 
ity in Fig. H is symmetric with respect to the x and y 
axes. In such a case, the wave functions can be divided 
into four symmetry classes 



(27) 
(28) 



with the parities £ € { — , +} and £ e { — , +}. The normal 
derivative obeys the same symmetry relations. 









/ r ^\ 




\ / x 







FIG. 2: Symmetric cavity. 

For systems with symmetries the BIEs can be reduced 
to a fundamental domain if a modified Green's function 
is used. This decreases the numerical effort considerably. 
Let us restrict our discussion to the case in Eqs. ((2^) and 
( p8| ) ; other symmetries can be treated in a similar way. 
The BIEs ( |l2| ) reduce to integrals along the boundaries 
restricted to the quadrant x, y > if the Green's function 
G(r,r') is replaced by 

G(r, r') + CG(n, r') + C£G(r 2 , r') + £G(r 3) r') (29) 

with r = (x,y). r 1 = (-x,y), r 2 = (-x,-y), r 3 

—y); see Fig. M. The derivative <9„G(s,r') is modified 
in the same way with the normal unit vector j/ changing 
as r. 



The scattering problem as formulated in Sec. II B does 
not allow the symmetry reduction because the incom- 
ing plane wave in general destroys the symmetry; <j>i n 
and ipi n in Eq. jl7| ) do not fulfil the conditions ( p7j ) and 
( p8| ) . There are certain incoming directions which do not 
spoil the symmetry, but using only these special direc- 
tions is dangerous because possibly not all resonances 
are excited. A better approach is to consider a different 
physical situation illustrated in Fig. |^. Four plane waves 
are superimposed to a symmetric incoming wave 

ipia = cxp (ikr)+Cexp (ikir)+CCexp (ik 2 r)+^exp (ik 3 r) 

(30) 

where k = {k x ,k y ), ki = (-k x ,k y ), k 2 = (-k x ,-k y ), 
k3 = (k x ,—k y ). With this incoming wave, the scatter- 
ing problem can be symmetry reduced. A more general 
formulation for an arbitrary symmetry can be found in 
Ref. [|33l . 



The most convenient numerical strategy for solving 
BIEs as in Eqs. (|§) and © is the BEM. The boundary 
is discretized by dividing it into small boundary elements. 
Along such an element, the wave function and its normal 
derivative are considered as being constant (for linear, 
quadratic, and cubic variations see, e.g., Refs. |2^, p6[). 
Equation ( |l3| ) is therefore approximated by a sum of Nj 
terms 



$3(B<J& + C«Vj) = 



(31) 



i=i 



where Bu = J t ds B(si,s), C« = ^dsC{si,s), 4>i = <f>(si), 
ipi = ip(si), and J, denotes the integration over a bound- 
ary element with midpoint s/. The entire set of BIEs 
is approximated by an equation as in Eq. (|l4|), but 
for which Bj and Cj are Nj x N matrices, M is a 
2N x 27V (non-Hermitian complex) matrix, ip and ip are 

A^-component vectors with 2N — Ylj=i^j- Note that 
each boundary element belongs to two different regions. 
In the same way the scattering problem is approximated 
by an equation as in Eq. (17) with Mo being a 2N x 27V 
matrix, <f>- in and ip- m being TV-component vectors. 

In the literature several levels of approximation are 
used for the matrix elements Bu and Ca. The crudest 
approximation is to evaluate such an integral only at the 
corresponding midpoint s/ . While this is sufficient for the 
calculation of bounded states in quantum billiards |p8| , 
in our case the small imaginary parts of k require a more 
accurate treatment. We therefore do perform the numer- 
ical integration of the matrix elements Bu and Ca , using 
standard integration routines like, for example, Gaussian 
quadratures [|34| . The number of interior points in the 
range of integration should be chosen large if the bound- 
ary elements s, and Si are close to each other and small 
if they are far away. Moreover, our experience is that the 
results are considerably more accurate if the boundary 
elements are not approximated by straight lines but, in- 
stead, the exact shape of the boundary elements is used 
for all interior points in the range of integration. 

Due to the almost singular behaviour of the integral 
kernels at r' = r, the diagonal elements Cu and Bu re- 
quire special care. Inserting the limiting cases for small 
boundary-element length As/ in Eqs. ( |l0| ) and ( |i"l| ) into 
Eq. (0) gives 



Ca 



-1 



2n 



(32) 



where k; is the curvature at point s/. To approximate 
Bu accurately, more higher order terms than in Eq. (|^) 
are needed: 



/i) / s 2i, z 2i 
#o z ~-ln- + l + -7, 

TT Z TT 



(33) 



where 7 = 0.577215 ... is Euler's constant. Integration 
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yields 



As, 



1 - In 



n-jkAsi 7r 

3 ; ' +1- 



•7] 



(34) 



Treatment of corners 



Dielectric corners are numerically difficult to treat be- 
cause certain components of the electric field can be in- 
finite at the corner (see the discussion in Ref. |i5| in the 
context of dielectric waveguides). In the BEM, a corner 
leads to a second problem. The integral kernel of C« 
has a singularity caused by a diverging curvature K] see 
Eq. (11). To circumvent these difficulties, we smooth the 
boundary as sketched in Fig. 0. The curvature k and the 
electric field are then everywhere bounded. 

The minimum value of the radius of curvature, p = 
1/k, along such a rounded corner should be much larger 
than the typical distance between discretization points, 
so that the boundary is locally smooth. However, in or- 
der to ensure that the rounding does not influence the 
result, p should be much smaller than the wavelength A. 
Clearly, these requirements can be met most efficiently by 
using a nonuniform discretization with a relatively large 
density of discretization points at corners as illustrated 
in Fig. ||. Since the results do not depend on the partic- 
ular selected rounding and discretization we do not give 
explicit formulae. 




FIG. 3: Rounded corner. The number of discretization points 
(circles) is enhanced at the corner. 



B. Finding and computing resonances 



The scattering problem as discussed in Sec. II B pro 



vides us with first approximations to the wave numbers 
of the resonances. Let us fix <p to an appropriate value 
and plot the total cross section in Eq. p0[ ) as function of 
k in the range of interest. Resonances can be identified 
as peaks. The peak position a and the width 7 deter- 
mine the resonant wave number as fc ros w fej = a — 27/2. 
It might be difficult to resolve numerically very broad 
and very narrow peaks, because they are hidden either in 
the background or between two consecutive grid points. 
For microlaser operation, however, these two extreme 
cases are not relevant. Too short-lived resonances (broad 
peaks) fail to provide a sufficient long lifetime for the 
light to accumulate the gain required to overcome the 



lasing threshold, whereas too long-lived resonances (nar- 
row peaks) do not supply enough output power. 

The spurious solutions of the interior Dirichlet prob- 
lem occasionally appear in the scattering spectrum as 
extremely narrow peaks. The reason is that numerical 
inaccuracies broaden the (5-peaks to peaks of finite width. 
However, choosing a sufficiently fine boundary discretiza- 
tion and/or an appropriate, not too fine discretization in 
k reduces the probability of observing them. Moreover, 
they can be removed with a simple trick: use k with a 
small negative imaginary part in Eq. (|2C|). 

The discretized version of Eq. (|lj) has a nontrivial so- 
lution only if det M(k ICS ) = 0. Using a first approxima- 
tion ki from the scattering problem as starting value, we 
find a much better approximation to k rcs in the complex 
plane with the help of Newton's method 



h+i = h- 



g'(ki) 



(35) 



with I = 1,2,... and g(k) = detM(fc). The derivative 
g'(k) = dg(k)/dk can be approximated by 



g'(k) 



g(k + A)-g(k) .g(k + iA)-g{k) 



2A 



2A 



(36) 



where A is a small real number. Equation ( p5| ) is re- 
peated iteratively until a chosen accuracy is achieved. 

Newton's method in Eq. ([35]) is very efficient close to 
an isolated resonance where detM oc k — k lcs . For q- 
fold degenerate resonances the determinant behaves like 
(k — k ICS ) q . The resulting problem of slow convergence 
can be eliminated by choosing g = (detM) 1 / 9 . 

A slightly different approach for finding resonances can 
be gained by rewriting Newton's method in Eq. ( |35| ) with 
the help of the matrix identity In det M = tr In M as 



k, 



1+1 



= k,~ 



(37) 



where tr denotes the trace of a matrix. The derivative 
M'(k) can be calculated as in Eq. (|3^). It turns out that 
the numerical algorithm corresponding to Eq. ( |37| ) is a 
bit faster than the original Newton's method in Eq. (|35|). 

Having found a particular wave number fc rcs , the vector 
components <pi and ipi are given by the null eigenvector of 
the square matrix M(k Ies ) . This eigenvector can be found 
with, for instance, singular value decomposition |3^| . The 
wave function in each domain Qj is then constructed by 
discretizing Eq. (^]) 

V>(r') =J2^ fdsd v G{ S) T';k ies ) (38) 

~^2<t>l / dsG{s,r';k ICS ) , 
I ^ l 

where I runs over all boundary elements of . 

How fine must be the discretization of the boundary 
in order to obtain a good approximation of a resonance 
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at /c res ? The local wavelength A = 27r/nRe(fc ros ) is the 
smallest scale on which the wave function and its deriva- 
tive may vary. Hence, the minimum number of bound- 
ary elements along each wavelength, b — A/ As, should 
be larger or equal than at least 4; As is the maximum 
value of all lengths As*. We have verified the BEM us- 
ing different values of b. Taking b = 16, we find good 
agreement with the separation-of- variables solution of the 
circular cavity (see e.g. Ref. |2(J) and to results of the 
wave-matching method obtained for the quadrupolar cav- 
ity (I). Only for extremely long-lived resonances larger 
b are necessary to determine the very small imaginary 
parts of k accurately (recall that this is important for 
distinguishing spurious solutions from real resonances). 
However, as already explained, extremely long-lived res- 
onances are not relevant for microlaser applications and, 
moreover, they occur only in circular or slightly deformed 
circular cavities for which the wave-matching method is 
more suitable anyway. 



IV. EXAMPLE: TWO COUPLED 
HEXAGONAL-SHAPED CAVITIES 

Vietze et al. have experimentally realized hexagonal- 
shaped microlasers by putting laser active dyes into 
molecular sieves made of AIPO4 — 5 H|. Numerical 
simulations on rounded hexagons based on the wave- 
matching method have shown convergence problems at 
corners ^6). The following example is relevant for 
future experiments and demonstrates that the BEM can 
handle arbitrarily sharp corners and, moreover, coupled 
resonators. Near-field-coupling of resonators is interest- 
ing, because it may improve the optical properties of the 
resonators, as e.g. the far-field directionality. 

Figure ^ illustrates the configuration: two hexag- 
onal cavities with sidelength R are displaced by the 
vector(1.8i?, 0.5R). According to the experiments in 
Ref. jlj, |l5| |, the polarization is of TM type, the index 
of refraction is n — 1.466 inside the cavities and n = 1 
outside; R ranges from 4/itm to 10/Ltm, the wavelength 
A from 600nm to 800nm depending on the dye. Since 
only the ratio between R and A is relevant, we use in the 
following the dimensionless wave number kR. We focus 
on a /ci?-interval from 20 to 25 within the experimental 
spectral interval. A total of 2N = 3200 discretization 
points is then sufficient. We slightly smooth the cor- 
ners as discussed in Sec. [II A such that p/X pa 0.11 and 
p/As ps 11.2. 

Figure || shows the total cross section a for plane-wave 
scattering with incidence angle (f> — 15° computed from 
Eq. (|2^). The dominant structure is a series of equidis- 
tant peaks of roughly Lorentzian shape. At kR pa 23.25 
we identify a spurious solution of the interior Dirichlet 
problem. The fact that it is the only one visible in the 
chosen range of wave numbers confirms that the spurious 
solutions are not a problem. 

The peak at kR pa 22.95 in Fig. || has roughly the 






x/R 



FIG. 4: Two hexagonal cavities. The incoming plane wave with 
wave vector k is incidence at 15° to the horizontal side faces. 




FIG. 5: Calculated total cross section a/R vs. kR for two coupled 
hexagonal resonators. The plane wave is incidence at 15° to the 
horizontal side faces; cf. Fig. W 



width 0.196, so we use k x R = 22.95 - i0.098 as initial 
guess for Newton's method in Eq. (37). The more pre- 
cise location of the resonance is found to be k Tes R ps 
22.94444 - i0.09696. The near-field intensity pattern in 



Fig. ^ and the far-field emission pattern in Fig. ^ are 
computed with the help of Eq. ([58). A detailed account 
of the structure of this kind of resonances and its impli- 
cation on the properties of the microlasers will be given 
in a future publication. 



V. SUMMARY 

We have introduced a boundary element method 
(BEM) to compute TM and TE polarized resonances 
with intermediate lifetimes in dielectric cavities. We have 
discussed spurious solutions, the treatment of cavities 
with symmetries and cavities with corners. Numerical 
results are shown for an example of two coupled hexago- 
nal cavities. 

If compared to finite-difference methods and related 
methods the BEM is very efficient since the wave function 
and its derivative are only evaluated at the boundaries 
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of the cavities. It is in general less efficient than the 
wave-matching method but in contrast to the latter it 
can be applied to complex geometries, such as cavities 
with corners and coupled cavities. 




FIG. 6: Calculated near-field intensity pattern |i/>(r)| 2 of the res- 
onance with fcrcs-R ~ 22.94444 - iO. 09696. Intensity is higher for 
light regions and vanishes in the black regions. 



The BEM is especially suitable for computing phase 
space representations of wave functions such as the 
Husimi function which also only requires the wave func- 
tion and its normal derivative on the domain bound- 
aries p7fl . 
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FIG. 7: Far-field emission pattern, \ip(r,d)\ 2 with large r, of the 
resonance in Fig. tl 



I would like to thank M. Hentschel, S. W. Kim, J. 
Nockel, F. Laeri and A. Backer for discussions. The work 
was supported by the Volkswagen foundation (project 
"Molekularsieblaser-Konglomerate im Infraroten" ) . 



[1] J. U. Nockel and A. D. Stone, Nature 385, 45 (1997). 

[2] C. Gmachl et al., Science 280, 1556 (1998). 

[3] A. Mekis et al, Phys. Rev. Lett. 75, 2682 (1995). 

[4] S. Chang, R. K. Chang, A. D. Stone, and J. U. Nockel, 

J. Opt. Soc. Am. B 17, 1828 (2000). 
[5] S. Lacey and H. Wang, Opt. Lett. 26, 1943 (2001). 
[6] J. U. Nockel and A. D. Stone, in Optical Processes in Mir- 

cocavities, Vol. 3 of Advanced Series in applied Physics, 

edited by R. Chang and A. Campillo (World Scientific, 

Singapore, 1995). 
[7] O. A. Starykh, P. R. J. Jacquod, E. E. Narimanov, and 

A. D. Stone, Phys. Rev. E 62, 2078 (2000). 
[8] S. Gianordoli et al., IEEE J. Quantum Electronics 36, 



458 (2000). 

[9] K. Shhna, R. Omori, and A. Suzuki, Opt. Lett. 26, 795 

(2001) . 

[10] S. B. Lee et al, Phys. Rev. Lett. 88, 033903 (2002). 
[11] N. B. Rex et al, Phys. Rev. Lett. 88, 094102 (2002). 
[12] M. Hentschel and K. Richter, Phys. Rev. E 66, 056207 

(2002) . 

[13] A. W. Poon, F. Courvoisier, and R. K. Chang, Opt. Lett. 

26, 632 (2001). 
[14] U. Vietze et al, Phys. Rev. Lett. 81, 4628 (1998). 
[15] I. Braun et al, Appl. Phys. B: Lasers Opt. 70, 335 (2000). 
[16] J. D. Jackson, Klassische Elektrodynamik (Walter de 

Gruyter, Berlin, New York, 1983). 



9 



[17] R. H. Landau, Quantum Mechanics II, 2 ed. (John Wiley 

& Sons, New York, 1996). 
[18] G. Gamow, Z. Phys. 51, 204 (1928). 
[19] P. L. Kapur and R. Peierls, Proc. Roy. Soc. Lond. A 166, 

277 (1938). 

[20] P. W. Barber and S. C. Hill, Light scattering by par- 
ticles: computational methods (World Scientific, Singa- 
pore, 1990). 

[21] P. M. van den Berg and J. T. Fokkema, IEEE Trans. 
Antennas Propag. 27, 577 (1979). 

[22] M. Lohmeyer, Opt. Quantum Electron. 34, 541 (2002). 

[23] G. Chen and J. Zhou, Boundary element methods (Aca- 
demic Press, San Diego, 1992). 

[24] M. Kitahara, Boundary integral equation methods in 
eigenvalue problems of elastodynamics and thin plates 
(Elsevier, Amsterdam, 1985). 

[25] Boundary element methods in acoustics, edited by R. D. 
Ciskowski and C. A. Brebbia (Computational Mechanics 
Publications and Elsevier Applied Science, Southampton 
Boston, 1991). 

[26] P. K. Banerjee, The boundary element methods in engi- 
neering (McGraw-Hill, London, 1994). 
[27] I. Kosztin and K. Schulten, Int. J. mod. Physics C 8, 293 



(1997). 

[28] A. Backer, e-print nlin.CD/0204061 (2002). 

[29] A. J. Burton and G. F. Miller, Proc. R. Soc. Lond. A 

323, 201 (1971). 
[30] P. Gaspard and S. A. Rice, J. Chem. Phys. 90, 2255 

(1989). 

[31] P. A. Knipp and T. L. Reinecke, Phys. Rev. B 54, 1880 
(1996). 

[32] I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, 
Series, and Products (Academic Press, New York, 1965). 

[33] J. U. Nockel and A. D. Stone, Phys. Rev. B 50, 17415 
(1994). 

[34] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. 
Vetterling, Numerical Recipes in C. The Art of Scientific 
Computing. (Cambridge University Press, Cambridge, 
1988). 

[35] G. R. Hadley, J. Lightwave Technol. 20, 1219 (2002). 
[36] J. U. Nockel, private communication (2002) (unpub- 
lished) . 

[37] M. Hentschel, H. Schomerus, and R. Schubert, e-print 
arXiv:physics/0208006 (2002). 



